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I.  \  ELLIPTIC  SOLVER  LIBRARY 

All  two  dimensional  Helmholtz  solvers  for  rectangles  at  NORDA 
have  been  placed  In  one  library  with  standardized  naming  and  argument 
list  conventions.  Each  solver  has  full  internal  documentation,  test 
software  for  accuracy,  timing  and  error  conditions  (such  as  incorrect 
dimensions,  etc.).  The  Internal  documentation  for  three  representative 
solvers  (SMNUUH,  SMNSUH  and  SMMSSH)  Is  reproduced. In  Appendix  A. 

To  Invoke  a  solver,  two  routines  must  be  called.  The  first 
(e.g.,  CMNUUH)  tests  the  correctness  of  the  arguments  and  produces 
constants  for  use  by  the  solver  proper;  It  need  be  called  only  once  for 
each  different  Helmholtz  operation  (l.e.,  it  does  not  depend  on  the 
right  hand  side  of  the  equation).  The  second  (e.g.,  SMNUUH)  when  called 
with  the  precalculated  constants  and  a  right  hand  side  array  actually 
supplies  the  solution  to  Helmholtz's  equation.  It  must  obviously  be 
called  once  for  each  solution  required.^  These  rectangular  solvers 
generally  must  have  the  same  boundary  condition  type  on  the  facing 
parallel  sides  and  uniform  grid  versions  are  generally  faster  than 
stretched  grid  solvers  so  the  following  convention  has  been  adopted.  If 
the  solver  Is  called  ABCDEF  then: 

A  -  Indicates  the  solver  routine  type.  It  can  be 
C  -  for  the  constant  setup  routine, 

S  -  for  the  acutal  solver  routine 

B  -  Indicates  the  boundary  condition  type  on  the  left  and  right 
sides;  It  can  be 
P  -  for  periodic, 

D  -  for  Dlrlchlet, 

N  -  for  (staggered  grid)  Neumann, 

M  -  for  mixed  (Dlrlchlet  and/or  Neumann). 


C  -  indicates  the  boundary  condition  type  on  the  top  and  bottom 
boundaries.  It  can  be  P,  D,  N  or  M. 

0  -  Indicates  the  grid  type  In  the  x-dl recti  on;  it  can  be 
U  -  for  uniform  grid  (constant  delta-x), 

S  -  for  stretched  grid  (variable  delta-x). 

E  -  indicates  the  grid  type  In  the  y-dlrectlon.  It  can  be  U  or 

$. 

F  -  indicates  the  elliptic  operator  type;  It  can  be 
H  -  for  Helmholtz's  equations, 

S  -  for  separable  elliptic  equations, 

P  -  for  general  positive  definite  elliptic  equations, 

G  -  for  general  (non-symmetrlc)  elliptic  equations. 

The  solvers  currently  In  the  library  are: 

•  SMNUUH 

•  SMDUUH 

•  SMPUUH 

e  SNMSUH 

e  SMDSUH 

•  SMPSUH 

•  SMMSSH. 

By  adding  a  driver  routine  with  appropriate  transpose  operations, 
solvers  SPMUUH,  SNMUSH,  SDMUSH  and  SPMUSH  can  effectively  be  created 
from  this  set  (the  appropriate  code  Is  contained  In  the  test  program). 

The  uniform-uniform  solvers  use  the  FACR(l)  algorithm  and 
therefore  restrict  the  active  x  (or  first)  dimension  to  be  odd  and  the  y 
dimension  to  be  of  the  form  'p'  times  2  to  the  power  ' q '  plus  V,  where 
'p*  Is  usually  3,  4  or  5  (although  7,  9,  11,  etc.  are  possible),  ' q '  Is 


a  positive  integer,  and  V  is  -1,  0  or  1  depending  on  the  boundary 
conditions.  They  are  the  fastest  solvers  In  the  library  and  need  the 
least  storage  for  constants.  The  stretched-unlform  solvers  use  the 
FACR(O)  algorithm,  which  is  very  similar  to  FACR(l)  but  places  no 
restrictions  on  the  x  dimension.  They  are  only  slightly  slower  then  the 
uniform-uniform  solvers  but  require  more  storage  for  constants. 
Stretched-unlform  Helmholtz  solvers  on  a  rectangle  are  equivalent  to 
uniform-uniform  solvers  In  the  surface  of  a  sphere,  and  therefore  find 
use  in  spherical  coordinate  ocean  models.  The  stretched-stretched 
solver  uses  matrix  eigenvalue  deposition  which  has  no  restrictions  on 
the  rectangle  dimensions  but  takes  a  long  time  to  generate  the 

relatively  large  number  of  required  constants.  More  importantly.  It  Is 
an  0  (n  cubed)  method  rather  than  the  0  (n  squared,  log  n)  performance 
of  the  FACR  algorithms,  and  therefore  becomes  relatively  slower  for 
larger  problems.  On  the  other  hand,  the  method  vectorizes  exceptionally 
well  and  Is  very  flexible.  It  Is  probably  the  best  method  available  for 
very  small  problems  (although  this  may  not  be  very  Important). 
Stretched-stretched  Helmholtz  solvers  can  be  applied  to  any  separable 

elliptic  equation  so  SWSSH  could  equally  have  been  called  SMMUUS  (or 
SMMSSS). 

Most  of  the  computation  within  the  solvers  Is  performed  by  two 
auxiliary  packages  which  might  be  useful  In  their  own  right  In  some 

applications.  The  first  Is  a  real  fast  Fourier  transform  (i.e.,  fast 
sine  and  cosine  transform)  package  by  JAYCOR  consultant  D.  R.  Moore. 
This  standard  FORTRAN  code  Is  suitable  for  automatic  vectorlzatlon  or 
vector  machines  and  Is  probably  the  fastest  transportable  FFT  package 
available.  It  contains  all  the  transforms  needed  by  FACR  solvers. 

Including  the  non-standard  transform  required  for  staggered  grid  Neumann 
boundaries,  and  Is  fully  documented  with  Its  own  test  routines.  The 
second  package  contains  routines  for  solving  various  sets  of  constant  or 
variable  coefficient  trldi agonal  systems  of  linear  equations.  It  also 
Is  written  In  standard  FORTRAN  to  allow  automatic  vectorlzatlon  and 
contains  full  Internal  documentation. 
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II.  NON-SEPARABLE  ELLIPTIC  SOLVERS 

The  existing  stabllzed  marching  code  was  found  to  be  the 
fastest  for  the  non-separable  equations  encountered  In  rigid  lid  ocean 
models  with  bottom  topography.  But  this  method  has  several 
disadvantages: 

e  The  method  Is  very  sensitive  to  the  placement  of  the  'block1 
boundaries  which  artificially  divide  up  the  region. 

e  Storage  and  time  grow  with  the  cube  of  the  dimension  of  a 
side  of  the  rectangle.  It  Is  therefore  best  suited  to  small 
or  mid-size  (e.g.,  75  x  150)  problems. 

•  Double  precision  arithmetic  must  be  used  throughout.  This 
obviously  doubles  the  storage  requirement  and  In  many  cases 

also  doubles  the  CPU  time.  On  the  CYBER  205  the  highest 

vectorlzable  precision  Is  only  accurate  to  about  14 
significant  figures,  so  marching  methods  are  at  a  large 

disadvantage  on  this  machine. 

A  diagnostic  printout  capability  has  been  added  to  the  package  to 
simplify  the  choice  of  block  boundary  positions.  The  other 
disadvantages  are  Intrlnlslc  to  the  method,  but  note  that  despite  them 
the  method  Is  still  the  fastest  available  (at  least  for  mid-size 
probl ems). 

The  best  Iterative  method  was  found  to  be  'Black  Box'  Multi  grid 
[Dend|y,  1982].  This  method  was  originally  developed  for  elliptic 
problems  with  discontinuous  coefficients;  the  coefficients  In  ocean 
model  cases  are  smooth  but  the  discontinuous  capability  allows  non- 

rectangular  problems  to  be  solved  by  Imbedding  In  a  rectangle  (with 
discontinuities  at  the  region  boundary).  Always  solving  In  rectangles 
greatly  simplifies  the  task  of  vector! zatl on,  and  defining  the  region 
via  the  coefficients  removes  the  necessity  for  the  complicated  data 
structures  associated  with  other  non- rectangular  multi-grid  codes. 
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However,  the  cost  of  this  simplification  Is  that  a  solution  is 
calculated  over  'land'  where  It  Is  not  needed  (increasing  solve  time), 
and  additional  storage  Is  required  because  all  course  grid  elliptic 
operators  are  9-polnt  even  though  the  fine  grid  operator  Is  5-polnt. 

The  total  storage  required  for  constants  is  7.4  meshes  and  one 
'Iteration'  takes  about  as  long  as  one  application  of  the  very  fast 
FACR(O)  solver  (which  only  works  on  rectangular  Helmholtz  problems)  and 
reduces  the  error  by  a  factor  of  10  to  50.  Only  three  or  four  such 
Iterations  are  required  when  used  in  ocean  models  so  the  method  Is  very 
fast.  Marching  methods  are  faster  for  small  problems  but  the  multi  grid 
code  has  the  advantage  that  both  storage  and  solve  time  Increase  only 
linearly  with  mesh  size.  Thus,  Black  Box  multi  grid  Is  the  method  of 
choice  for  large  problems,  and  may  even  be  competitive  for  solving 
Helmholtz  equation  in  non-rectangular  regions  (despite  the  existence  of 
very  fast  direct  Helmholtz  solvers). 
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MODEL  COMPARISON 


A  regional  eddy -mean  energetics  package  has  been  developed  for 
the  quasi -geostroph 1c  ocean  model.  This  was  tested  on  a  course  grid  (25 
x  25)  version  of  a  single  gyre  experiment  reported  on  by  Holland 
[Holland,  1978].  Good  agreement  was  found  and  sample  output  from  the 
package  Is  reproduced  In  Appendix  8. 

The  series  of  wind  driven  mesoscale  eddy  experiments  for  the 
model  comparison  project  were  completed  and  the  results  presented  at  the 
1982  Fall  AGU  Meeting.  The  basic  finding  was  that  In  a  classic  double 
gyre  experiment  the  mid-latitude  jet  develops  much  more  slowly  In  the 
quasl-geostrophlc  results.  There  were  also  significant  differences  In 
the  edqy-mean  energetics. 
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APPENDIX  A 


Sample  Internal  Documentation  from  Rectangular  Helmholtz  Solver  Library 

•  CMNUUH  and  SMNUUH 

•  CMNSUH  and  SMNSUH 

•  CMMSSH  and  SMMSSH 


SUBROUTINE  CMNUUH(DAX,DAY,DAC,  NDX , NX , N Y ,  Q ,  LQ ,  U) 

INTEGER  NDX, NX, NY,  LQ 

DOUBLE  PRECISION  DAX ( 3, 1 > , DAY < 3 , 1 > , DAC 
DOUBLE  PRECISION  U 
DIMENSION  Q(LO>,  W<NY,5> 

C 

CAAAAAAAAAAAAAAAAAAAAAAAAAA  AAAAAAAAAA  A  A  A  A  A  AAA  AAA  A  A  AAA  A  AAAAA AA  AAA  AAAAAA  A  A 
C AAAAAAAAAA A A AAA AAA AAA AAA Ax A AAA A A AAAAA A AAAAAAAA AAAAAAAAA A AAAA A AAA A A A A * A  A 


CAA 

AA 

CAA 

INTRODUCTION: 

AA 

CAA 

AA 

CAA 

COMPUTES 

THE  CONSTANTS  REQUIRED  BY  SMNUUH  TO  SOLVE  THE 

AA 

CAA 

FOLLOWING  <  MIXED-NEUMANN  )  HELMHOLTZ ' S  EQUATION  ON 

AA 

CAA 

AN  NXANY 

MESH. 

AA 

CAA 

AA 

CAA 

FOR  I  * 

1,  ...»  NX  AND  3  =  1,  ...,  NY 

AA 

CAA 

AA 

CAA 

AN(  J)  A  H < I,  3+1)  + 

AA 

CAA 

AS (  J)  A  H ( I ,  3-l>  ♦ 

AA 

CAA 

AE< I  )  A  H(I+1,J  )  + 

AA 

CAA 

AW<  I  )  A  H  <  1-1 ,  J  >  + 

AA 

CAA 

AP< I, J)  A  H< I,  J  >  •  S( I,J> 

AA 

CAA 

AA 

CAA 

WHERE 

AA 

CAA 

AN  ( J  ) 

»  DAY (3,1)  3*1. ..NY-1, 

AA 

CAA 

0  3*  NY. 

AA 

CAA 

AS  ( J  ) 

*  0  3*  1. 

AA 

CAA 

DAY (3,1)  3*2  ...  NY , 

AA 

CAA 

AE(  I) 

*  DAX (3,1)  1*1. ..NX-1, 

AA 

CAA 

0  I*  NX. 

AA 

CAA 

AW  (  I  > 

*  0  I*  1. 

AA 

CAA 

DAX (3,1)  1*2. ..NX, 

AA 

CAA 

AP  ( I , 

3)  *  DAC+DAX<2,1 )-< AN(3)+AS(3)*BE( I)+BW( I) > 

AA 

CAA 

1*1. ..NX,  3*1. ..NY. 

AA 

CAA 

BE(  I) 

-  DAX (3,1)  1*1. ..NX-1, 

AA 

CAA 

DAX < 1 , 1 )  I*  NX. 

AA 

CAA 

BW<  I) 

*  DAX (1,1)  I*  1. 

AA 

CAA 

DAX (3,1)  1*2. ..NX, 

AA 

CAA 

AA 

CAA 

USUALLY 

DAX ( 2 , 1 ) *0 . DO  . 

AA 

CAA 

AA 

CAA 

FOR  DIRICHLET-NEUMANN  PROBLEMS  DAX ( 1 , 1 > «DAX ( 3 , 1 )  , 

AA 

CAA 

FOR  NEUMANN-NEUMANN  PROBLEMS  DAX ( I , 1 ) =0 .DO  . 

AA 

CAA 

AA 

CAA 

FOR  COMPATABILITY  WITH  OTHER  SOLVERS  SET  DAY< 1 , I >»0.D0 

AA 

CAA 

AND  DAY ( 2 , 1 >  *0 . DO  . 

AA 

CAA 

AA 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

CAA 

AA 

CAA 

ARGUMENT 

LIST:  SEE  OTHER  COMMENTS  FOR  DETAILS 

AA 

CAA 

AA 

CAA 

DAX, DAY, 

DAC  -  HELMHOLTZ  COEFFICIENTS 

AA 

CAA 

AA 

CAA 

NDX 

-  1ST  ARRAY  DIMENSION  OF  S  AND  H  IN  THE  SOLVER 

AA 

CAA 

AA 

CAA 

NX, NY 

-  DIMENSIONS  OF  MESH  (NX  ODD  AND  NY-PA2AAL) 

AA 

CAA 

AA 

CAA 

Q(LQ> 

-  ARRAY  FOR  REQUIRED  CONSTANTS 

AA 

CAA 

AA 

CAA  LQ  -  SIZE  OF  Q,  MUST  BE  AT  LEAST  LQ1+LQ2+5,  WHERE:  AA 
CAA  LQ1  «  NY  +  10ALQ3AA2  ♦  1  AA 
CAA  LQ2  ■  NYA( <NX+7>/4>  *  <NY*3>/2  +2  AA 
CAA  LQ3  -  0  IF  P  *  4,5,  OR  6;  AA 


CAA  <P-l>/2  IF  P  =  7  OR  GREATER  AA 
CAA  A* 
CAA  W  (  N  Y  ,  5  >  -  DOUBLE  PRECISION  WORKSPACE  ARRAY  AA 
CAA  AA 
CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


DOUBLE  PRECISION  WORKSPACE  ARRAY 


SPECIAL  FEATURES  OR  NON-ANSI  USAGE 

AN  ERROR  MESSAGE  WILL  BE  OUTPUT  TO  CHANNEL  & ,  AND  THE  PROGRAM 
TERMINATED,  IF  ANY  OF  THE  FOLLOWING  CONDITIONS  ARE  DETECTED: 
NX  LESS  THAN  5  OR  EVEN 

NY  INCOMPATABLE  WITH  THE  FOURIER  PACKAGE 

LQ  TOO  SMALL 

DAX,DAY,DAC  NOT  CONSISTANT  WITH  HELMHOLTZ ' S  EQUATION 

THE  WORKSPACE  ARRAY,  W,  NEED  NOT  BE  DISTINCT  FROM  THE 
WORKSPACE  SUPPLIED  TO  THE  SOLVER.  HOWEVER  THE  LATTER  MAY 
CONSIST  OF  SINGLE  PRECISION  ARRAYS. 

FOR  MORE  DETAILED  DOCUMENTATION  -  SEE  SMNUUH. 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA A AAA AAAAAAAAAAAAAAAAAAAAAAAAAAA 
CA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
C 


C  END  OF  CMNUUH. 

END 
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SUBROUTINE  SMNUUH ( S , H , NDX , NX , NY  ,  Q.LQ,  U1,W2> 

INTEGER  NDX, NX, NY,  LQ 
DIMENSION  S(NDX,NY> ,H(NDX,NY) 

DIMENSION  Q(LQ> 

DIMENSION  Ul(l),  W2(l> 

C . DIMENSION  U1 (NX+1,NY> ,W2(NX,NY> 

C 

CAAA ************************ AAA A  AAA  *************  A  AAA*  ******************  A 
CAAAAAAAAAAAAAAAAAAAAAAAAAAA**********^ I********************************* 


CA* 

AA 

CAA 

DATE:  SEPTEMBER  1978 

AA 

CAA 

AA 

CAA 

AUTHOR:  D.R.  MOORE,  A.J.  UALLCRAFT 

AA 

CAA 

AA 

CAA 

ARCHIVIST:  A.J.  UALLCRAFT 

AA 

CAA 

AA 

CAA 

ADDRESS:  NORDA 

AA 

CAA 

CODE  322 

AA 

CAA 

NSTL  STATION,  MS  39529 

AA 

CAA 

AA 

CAA 

PHONE:  601-688-4835  (COMMERCIAL) 

AA 

CAA 

494-4835  (FTS) 

AA 

CAA 

485-4835  (AUTOVON) 

AA 

CAA 

AA 

CAA 

DATE  OF  CURRENT  SEQUENCE:  DECEMBER  1982 

AA 

CAA 

AA 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

CAA 

AA 

CAA 

INTRODUCTION: 

AA 

CAA 

AA 

CAA 

SOLVES  EINITE  DIFFERENCE  HELMHOLTZ'S  EQUATIONS: 

AA 

CAA 

OVER  A  RECTANGLE  IN  CARTESIAN  COORDINATES, 

AA 

CAA 

ON  A  UNIFORM-UNIFORM  MESH, 

AA 

CAA 

WITH  (ZERO)  NEUMANN  BOUNDARY  CONDITIONS  1/2  GRID  DISTANCE 

AA 

CAA 

OUTSIDE  THL  MESH  AT  J»1  AND  J»NY, 

AA 

CAA 

AND  EITHER  (ZERO)  NEUMANN  BOUNDARY  CONDITIONS  1/2 

GRID 

AA 

CAA 

DISTANCE  OUTSIDE  THE  MESH  OR  (ZERO)  DIRICHLET  BOUNDARY 

AA 

CAA 

CONDITIONS  1  GRID  DISTANCE  OUTSIDE  THE  MESH  AT  1= 

1  AND  I=NX. 

AA 

CAA 

AA 

CAA 

USES  THE  (UNSTABALISEID  FACR(l)  METHOD. 

AA 

CAA 

AA 

CAA 

BOTH  THE  INPUT  ARRAY  S(NDX,NY>  AND  THE  OUTPUT  ARRAY 

H ( NDX , NY ) 

AA 

CAA 

ARE  ASSUMED  TO  CONTAIN  ACTIVE  SUBARRAYS  FROM  1  TO  NX  AND 

AA 

CAA 

1  TO  NY.  THE  NON-ACTIVE  SUBARRAY,  NX*1  TO  NDX  AND 

1  TO  NY, 

AA 

CAA 

OF  S  CAN  BE  ARBITRARY  ON  INPUT  AND  THAI  OF  H  IS 

RETURNED 

AA 

CAA 

INTACT. 

AA 

CAA 

AA 

CAA 

IN  THE  CALLING  SEQUENCE  THE  INPUT  AND  OUTPUT  ARRAYS 

MAY  BE  THE 

AA 

CAA 

SAME  ARRAY.  IF  THEY  ARE  NOT  THE  SAME,  THE  INPUT  ARRAY  IS 

AA 

CAA 

RETURNED  INTACT  BY  THE  SOLVER. 

AA 

CAA 

AA 

CAA 

IT  IS  ASSUMED  THAT  NECESSARY  CONSTANTS  HAVE  BEEN  PRECOMPUTED 

AA 

CAA 

EXTERNALLY  AND  ARE  STORED  IN  THE  ARRAY  Q ( LQ ) ,  WHICH 

AA 

CAA 

IMPLICITLY  DEFINES  THE  HELMHOLTZ  COEFFICIENTS  (DAX, 

DAY,  DAC) 

AA 

CAA 

TO  THE  SOLVER. 

AA 

CAA 

AA 

CAA 

AA 

CAA 

FOR  I  ■  1,  ...,  NX  AND  J  -  1,  ...,  NY  THE  SOLUTION  WILL 

AA 

CAA 

SATISFY: 

AA 

CAA 

AA 

CAA 

AN (  J)  A  H( I,  J+l >  ♦ 

AA 

CAA 

AS (  J)  *  H( I,  J-l >  + 

AA 

CAA 

AE( I  )  A  H ( 1*1, J  >  + 

AA 
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CAA 

AW (  I  )  A 

H  < I- 1 , J  )  + 

AA 

CAA 

AP  < I , J )  A 

H(  I,  J  >  * 

S( I, J) 

A* 

CAA 

A  A 

CAA 

WHERE 

AA 

CAA 

AN  <  J ) 

*  DAY (3,1) 

J*1 . . .NY-1 , 

AA 

CAA 

0 

J*  NY. 

AA 

CAA 

AS  ( J  ) 

*  0 

J*  1. 

AA 

CAA 

DAY (3, 1 ) 

J*2. . .NY, 

AA 

CAA 

AE(  I) 

*  DAX (3,1) 

1*1. ..NX-1, 

AA 

CAA 

0 

I*  NX. 

AA 

CAA 

AW(  I) 

-  0 

I*  1. 

AA 

CAA 

DAX (3,1) 

1*2. ..NX, 

AA 

CAA 

AP< I, J> 

*  DAC+DAX ( 2 

,1)-<AN( J)+AS( J)+BE( I)+BW( I)  ) 

A  A 

CAA 

1*1. ..NX,  J*1 . . . NY  . 

AA 

CAA 

BE<  I) 

*  DAX (3,1) 

1*1. ..NX-1, 

AA 

CAA 

DAX (1,1) 

I*  NX. 

AA 

CAA 

BW<  I) 

*  DAX (1,1) 

I*  1. 

AA 

CAA 

DAX (3,1) 

1*2. . .NX, 

AA 

CAA 

AA 

CAAAAAAAAAAAAAA AA AAA AAA AAA AAAAA AAA A A A A  A A A* AAA AAAAAAAAAAAAAAAA AAAAAAAAAAA 
CAA  AA 

CAA  LIBRARIES  REQUIRED:  AA 

CAA  AA 

CAA  FAST  FOURIER  TRANSFORM  LIBRARY  BY  D.R.  MOORE  AA 

CAA  TRIDIAGONAL  LINEAR  EQUATION  SOLVER  LIBRARY  BY  A.J.  UALLCRAFT  AA 

CAA  AA 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA A A AAAAAAAAAAAAA AAA A AAA AAA A A 
CAA  AA 

CAA  ARGUMENT  LIST:  SEE  OTHER  COMMENTS  FOR  DETAILS  AA 

CAA  AA 

CAA  S ( NDX  f  NY )  -  INPUT  ARRAY  AA 

CAA  AA 

CAA  H<NDX,NY)  -  OUTPUT  ARRAY  A A 

CAA  AA 

CAA  NDX, NY  -  SIZES  OF  S  AND  H  AA 

-  ACTIVE  SUBDIMENSION  OF  S  AND  H 

-  CONSTANT  ARRAY,  PRECOMPUTED  EXTERNALLY 

-  DIMENSION  OF  Q 


CAA 

CAA  NX 
CAA 

CAA  Q(LQ) 
CAA 

CAA  LQ 


CAA  U1 < NX  +  1 , NY )  -  WORKSPACE  ARRAYS,  AA 
CAA  W2(NX,  NY)  RETURNED  AND  AVAILABLE  AFTER  EACH  CALL.  AA 
CAA  AA 
CAAAAAAAAAAAA AAAAAAA AA AA AAAAAAAAAAAAA AAA AA AAAAAAAAAAAAAAAAAA AAAAAAAAAA AA 
CAA  AA 
CAA  SPECIAL  FEATURES  OR  NON-ANSI  USAGE:  AA 
CAA  AA 
CAA  DOUBLE  PRECISION  LIBRARIES  USE  THE  IMPLICIT  STATEMENT  AA 
CAA  AA 
CAA  TI-ASC  ONLY:  USE  W-OPIION  WHEN  COMPILING  FOURIER  PACKAGES  AA 
CAA  AA 
CAA  NX  MUST  BE  ODD  AND  .GE.  5,  AA 
CAA  NY  MUST  BE  OF  THE  FORM  PA2AAL,  WITH  L  *  1,  ...»  7.  AA 
CAA  USUALLY  P*4,  5,  OR  &  -  AA 
CAA  HIGHER  VALUES  OF  P  ARE  POSSIBLE,  BUT  ARE  LESS  EFFICIENT.  AA 
CAA  LQ  MUST  BE  AT  LEAST  LQ1+LQ2+5  (SEE  CMNUUH).  AA 
CAA  AA 
CAA  PROBLEMS  WITH  NON-ZERO  NEUMANN  OR  DIRICHLET  BOUNDARY  CONDITIONS  AA 
CAA  CAN  BE  SOLVED  BY  SUITABLY  MODIFYING  THE  BOUNDARY  SOURCE  TERMS.  AA 
CAA  ( I . E .  MODIFY  S(I,J>  WHERE  1-1,  OR  I*NX,  OR  J*l,  OR  J*NY>  AA 
CAA  AA 
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I 


CAA  MUST  CALL  CMNUUH  FIRST  TO  SET  UP  CONSTANTS  REQUIRED  BY  SMNUUH.  A* 
CAA  IF  NONE  OF  (NX.NY.DAX.DAY.DAC)  CHANGE ,  THEN  SMNUUH  MAY  BE  Aa 
CAA  CALLED  REPEATEDLY.  IF  ANY  OF  THESE  CHANGE.  CMNUUH  MUST  BE  AA 
CAA  CALLED  TO  RECALCULATE  THE  REQUIRED  CONSTANTS.  IF  TWO  OR  MORE  AA 
CAA  EQUATIONS  ARE  SOLVED  REPEATEDLY.  ONE  AFTER  THE  OTHER.  THEN  ONLY  AA 
CAA  ONE  CALL  TO  CMNUUH  FOR  EACH  IS  REQUIRED  IF  THEY  ARE  ALLOCATED  AA 
CAA  DISTINCT  CONSTANT  ARRAYS,  Q.  AA 
CAA  AA 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

C 

C  WORK  IS  DONE  BY  SUBROUTINE  SMNUU1 . 

C 

MX*  <  NX  +  1 ) /2 
IW*MXANY*1 

CALL  SMNUU1 (S.H.NDX.NX.NY,  Q , LQ ,  U1 , U1 ( IU ) , U2, MX , U2 < IU > , MX-1 > 
RETURN 

C  END  OF  SMNUUH. 

END 


V 


SUBROUTINE  CMNSUH <  DAX , DAY , LAC ,  N  D  X  ,  N  X  ,  N  Y  ,  Q.LQ,  U) 

INTEGER  NDX,NX,NY,  LQ 

DOUBLE  PRECISION  DAX ( 3 , NX  )  , DAY (3,1), DAC 
DOUBLE  PRECISION  U 
DIMENSION  Q ( LQ  )  v  W(NY,5) 

C 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


CAA 

AA 

CAA 

INTRODUCTION: 

A  A 

CAA 

A  A 

CAA 

COMPUTES 

THE  CONSTANTS  REQUIRED  BY  SMNSUH  TO  SOLVE  THE 

AA 

CAA 

FOLLOWING  <  MIXED-NEUMANN  )  HELMHOLTZ'S  EQUATION  ON 

AA 

CAA 

AN  NXANY 

MESH. 

AA 

CAA 

AA 

CAA 

FOR  I  * 

1 ,  . . . ,  NX  AND  J  *  1 ,  . . . ,  NY 

AA 

CAA 

AA 

CAA 

AN  <  J)  A  H  < I,  J+l)  + 

AA 

CAA 

AS <  J)  A  H< I ,  J-l>  ♦ 

AA 

CAA 

AE ( I  )  A  H< 1+1 , J  )  + 

AA 

v 

CAA 

AW< I  >  A  H< 1-1, J  >  + 

AA 

"  * 

CAA 

AP< I, J)  A  H( I,  J  )  *  S< I,  J) 

AA 

CAA 

AA 

CAA 

WHERE 

AA 

CAA 

AN<  J) 

■  DAY (3,1)  J*l... NY-1, 

AA 

CAA 

0  J*  NY. 

AA 

CAA 

AS  ( J ) 

-  0  J*  1. 

AA 

CAA 

DAY  <  3 , 1 >  J*2...NY, 

AA 

•V 

CAA 

AE(  I) 

*  DAX (3,1)  1-1... NX-1, 

AA 

CAA 

0  I«  NX. 

AA 

CAA 

AW<  I) 

*  0  I*  1 . 

AA 

1-4 

CAA 

DAX (1,1)  1*2. ..NX, 

AA 

CAA 

AP  (  I , 

J)  *  DAC+DAX(2,1)-(AN(J)+AS(J)+BE(I>  +BW< I) > 

AA 

CAA 

1*1. ..NX,  J-1...NY. 

AA 

CAA 

BE  (  I) 

*  DAX(3, I)  1*1. ..NX. 

AA 

CAA 

BW(  I) 

*  DAX (1,1)  1*1. ..NX. 

AA 

CAA 

AA 

CAA 

USUALLY 

DAX ( 2 , I ) *0 . DO  . 

AA 

m 

CAA 

AA 

CAA 

FOR  DIRICHLET-NEUMANN  PROBLEMS  DAX < 1 , 1 > , DAX < 3 , NX >  ARE 

AA 

*- 

CAA 

SIMILAR  TO  OTHER  DAX  TERMS. 

A  A 

CAA 

FOR  NEUMANN-NEUMANN  PROBLEMS  DAX <1 , 1 > *DAX ( 3 , NX ) *0 . DO  . 

AA 

CAA 

AA 

CAA 

FOR  COMPATAB IL ITY  WITH  OTHER  SOLVERS  SET  DAY < 1 , 1 ) *0 . DO 

AA 

* 1 

CAA 

AND  DAY <  2 , 1 >*0 .DO  . 

AA 

CAA 

AA 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

CAA 

AA 

CAA 

ARGUMENT 

LIST:  SEE  OTHER  COMMENTS  FOR  DETAILS 

AA 

CAA 

AA 

CAA 

DAX, DAY , 

DAC  -  HELMHOLTZ  COEFFICIENTS 

AA 

CAA 

AA 

CAA 

NDX 

-  1ST  ARRAY  DIMENSION  OF  S  AND  H  IN  THE  SOLVER 

AA 

CAA 

AA 

CAA 

NX, NY 

-  DIMENSIONS  OF  MESH  (NY-PA2AAL  -  SEE  SMNSUH) 

AA 

CAA 

AA 

CAA 

Q  ( LQ ) 

-  ARRAY  FOR  REQUIRED  CONSTANTS 

AA 

CAA 

AA 

-* 

CAA 

LQ 

-  SIZE  OF  Q,  MUST  BE  AT  LEAST  LQ1+LQ2,  WHERE: 

AA 

CAA 

LQ1  ■  NY  *  10ALQ3AA2  ♦  1 

AA 

CAA 

LQ2  *  NXA ( NY  +  3  ) 

AA 

CAA 

LQ3  *0  IF  P  *  4,5,  OR  6 ; 

AA 

. 

CAA 

<P-l)/2  IF  P  *  7  OR  GREATER 

AA 

* 
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CAA  AA 
CAA  U ( NY , 5 )  -  DOUBLE  PRECISION  WORKSPACE  ARRAY  AA 
CAA  AA 
CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


CAA  AA 

CAA  SPECIAL  FEATURES  OR  NON-ANSI  USAGE  AA 


CAA  AA 
CAA  AN  ERROR  MESSAGE  WILL  BE  OUTPUT  TO  CHANNEL  6,  AND  THE  PROGRAM  AA 
CAA  TERMINATED,  IF  ANY  OF  THE  FOLLOWING  CONDITIONS  ARE  DETECTED:  AA 
CAA  NX  LESS  THAN  3  AA 
CAA  NY  INCOMPATABLE  WITH  THE  FOURIER  PACKAGE  AA 
CAA  LQ  TOO  SMALL  AA 
CAA  DAX,DAY,DAC  NOT  CONSISTANT  WITH  HELMHOLTZ'S  EQUATION  AA 
CAA  AA 


CAA  THE  WORKSPACE  ARRAY,  U,  NEED  NOT  BE  DISTINCT  FROM  THE  AA 
CAA  WORKSPACE  SUPPLIED  TO  THE  SOLVER .  HOWEVER  THE  LATTER  MAY  AA 
CAA  CONSIST  OF  SINGLE  PRECISION  ARRAYS.  AA 
CAA  AA 
CAA  FOR  MORE  DETAILED  DOCUMENTATION  -  SEE  SMNSUH.  AA 
CAA  AA 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAA 
C 


C 


END  OE  CMNSUH 
END 


SUBROUTINE  SHNSUH < S , H f NDX , NX , NY ,  Q.LQ,  W1,U2> 

INTEGER  NDX , NX , N Y ,  LQ 
DIMENSION  S(NDX,NY> ,H(NDX,NY > 

DIMENSION  Q( LQ ) 

DIMENSION  U1 (NDX, NY) ,W2(NDX,NY) 

C 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


CAA 

AA 

CAA 

DATE:  SEPTEMBER  1978 

AA 

CAA 

A  A 

CAA 

AUTHOR:  D.R.  MOORE,  A.J.  UALLCRAFT 

AA 

CAA 

A  A 

CAA 

ARCHIVIST:  A.J.  UALLCRAFT 

AA 

CAA 

AA 

CAA 

ADDRESS:  NORDA 

AA 

CAA 

CODE  322 

AA 

CAA 

NSTL  STATION,  MS  39529 

A  A 

CAA 

A  A 

CAA 

PHONE:  601-688-4835  (COMMERCIAL) 

AA 

CAA 

494-4835  (FTS) 

AA 

CAA 

485-4835  (AUTOVON) 

AA 

CAA 

AA 

CAA 

DATE  OF  CURRENT  SEQUENCE:  DECEMBER  1982 

AA 

CAA 

AA 

CAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

CAA 

AA 

CAA 

INTRODUCTION: 

AA 

CAA 

AA 

CAA 

SOLVES  FINITE  DIFFERENCE  HELMHOLTZ'S  EQUATIONS: 

AA 

CAA 

OVER  A  RECTANGLE  IN  CARTESIAN  COORDINATES, 

AA 

CAA 

ON  A  STRETCHED-UNIFORM  MESH, 

AA 

CAA 

WITH  (ZERO)  NEUMANN  BOUNDARY  CONDITIONS  1/2  GRID  DISTANCE 

AA 

CAA 

OUTSIDE  THE  MESH  AT  J*1  AND  J«NY, 

AA 

CAA 

AND  EITHER  (ZERO)  NEUMANN  BOUNDARY  CONDITIONS  1/2 

GRID 

AA 

CAA 

DISTANCE  OUTSIDE  THE  MESH  OR  (ZERO)  DIRICHLET  BOUNDARY 

AA 

CAA 

CONDITIONS  1  GRID  DISTANCE  OUTSIDE  THE  MESH  AT  I * 

1  AND  I*NX. 

AA 

CAA 

AA 

CAA 

OR  ANY  EQUIVALENT  SET  OF  EQUATIONS;  FOR  EXAMPLE: 

AA 

CAA 

HELMHOLTZ'S  EQUATION  IN  (SOME)  NON-CARTESIAN  COORDINATE 

AA 

CAA 

SYSTEMS,  OR  MORE  GENERAL  (SEPARABLE)  ELLIPTIC  EQUATIONS. 

AA 

CAA 

AA 

CAA 

USES  THE  FACR(O)  METHOD. 

AA 

CAA 

AA 

CAA 

BOTH  THE  INPUT  ARRAY  S(NDX,NY)  AND  THE  OUTPUT  ARRAY 

H(NDX.NY) 

AA 

CAA 

ARE  ASSUMED  TO  CONTAIN  ACTIVE  SUBARRAYS  FROM  1  TO  NX  AND 

AA 

CAA 

1  TO  NY.  THE  NON-ACTIVE  'SUBARRAY,  NX+1  TO  NDX  AND 

1  TO  NY, 

AA 

CAA 

OF  S  CAN  BE  ARBITRARY  ON  INPUT  AND  THAT  OF  H  IS 

RETURNED 

AA 

CAA 

INTACT. 

AA 

CAA 

AA 

CAA 

IN  THE  CALLING  SEQUENCE  THE  INPUT  AND  OUTPUT  ARRAYS 

MAY  BE  THE 

AA 

CAA 

SAME  ARRAY.  IF  THEY  ARE  NOT  THE  SAME,  THE  INPUT  ARRAY  IS 

AA 

CAA 

RETURNED  INTACT  BY  THE  SOLVER. 

AA 

CAA 

AA 

CAA 

IT  IS  ASSUMED  THAT  NECESSARY  CONSTANTS  HAVE  BEEN  PRECOMPUTED 

AA 

CAA 

EXTERNALLY  AND  ARE  STORED  IN  THE  ARRAY  Q ( LQ ) ,  WHICH 

AA 

CAA 

IMPLICITLY  DEFINES  THE  HELMHOLTZ  COEFFICIENTS  (DAX, 

DAY,  DAC) 

AA 

CAA 

TO  THE  SOLVER. 

AA 

CAA 

AA 

CAA 

AA 

CAA 

FOR  I  »  1,  ...,  NX  AND  J  «  1,  ...,  NY  THE  SOLUTION  WILL 

AA 

CAA 

SATISFY: 

AA 

CAA  AA 


CAA 

AN<  J)  A  H  <  I ,  0  +  l>  ♦ 

AA 

CAA 

AS (  J)  A  H  (  I ,  0-1)  + 

AA 

CAA 

AE< I  )  A  H ( 1+1,0  >  + 

kk 

CAA 

AW< I  )  A  H < 1-1 , J  )  + 

kk 

CAA 

AP< I, J)  A  H< I,  0  >  *  S< I, J) 

kk 

CAA 

kk 

CAA 

WHERE 

kk 

CAA 

AN<0) 

■  DAY (3,1)  0*1... NY-1, 

kk 

CAA 

0  0*  NY. 

kk 

CAA 

AS  ( 0  > 

■  0  0  =  1 . 

kk 

CAA 

DAY (3,1)  0*2 ...  NY , 

kk 

CAA 

AE<  I) 

*  DAX (3,1)  1*1. ..NX-1, 

kk 

CAA 

0  I*  NX. 

kk 

CAA 

AW<  I) 

*  0  I*  1. 

kk 

CAA 

DAX (1,1)  1*2. ..NX, 

kk 

CAA 

AP< I, J) 

*  DAC+DAX<2,1 >-< AN(0)+AS(0>+BE< I)+BW( I> ) 

kk 

CAA 

1*1. ..NX,  0*1. ..NY. 

kk 

CAA 

BE  ( I ) 

*  DAX ( 3 , I >  I-1...NX. 

kk 

CAA 

BW<  I) 

*  DAX (1,1)  1*1. ..NX. 

kk 

CAA 

kk 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

CAA 

kk 

CAA 

LIBRARIES  REQUIRED: 

kk 

CAA 

kk 

CAA 

EAST  FOURIER  TRANSFORM  LIBRARY  BY  D.R.  MOORE 

kk 

CAA 

TRIDIAGONAL 

LINEAR  EQUATION  SOLVER  LIBRARY  BY  A.O.  WALLCRAFT 

kk 

CAA 

kk 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAA AAA AAA AAAAAAAAAAAAAAAAAAAAA 

CAA 

kk 

CAA 

ARGUMENT  LIST:  SEE  OTHER  COMMENTS  FOR  DETAILS 

kk 

CAA 

kk 

CAA 

S(NDX,NY) 

-  INPUT  ARRAY 

kk 

CAA 

kk 

CAA 

H<NDX,NY) 

-  OUTPUT  ARRAY 

kk 

CAA 

kk 

CAA 

NDXf  NY 

-  SIZES  OF  S  AND  H 

kk 

CAA 

kk 

CAA 

NX 

-  ACTIVE  SUBDIMENSION  OF  S  AND  H 

kk 

CAA 

kk 

CAA 

Q(LQ> 

-  CONSTANT  ARRAY,  PRECOMPUTED  EXTERNALLY 

kk 

CAA 

kk 

CAA 

LQ 

-  DIMENSION  OF  Q 

kk 

CAA 

kk 

CAA 

W1 ( NDX,NY) 

-  WORKSPACE  ARRAYS, 

kk 

CAA 

W2 <  NDX ,  NY  > 

RETURNED  AND  AVAILABLE  AFTER  EACH  CALL. 

kk 

CAA 

kk 

C A AAA AAA AAA AAA AAA AAAAAA AAA AAAA AAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAA 

CAA 

kk 

CAA 

SPECIAL  FEATURES  OR  NON-ANSI  USAGE: 

kk 

CAA 

kk 

CAA 

DOUBLE  PRECISION  LIBRARIES  USE  THE  IMPLICIT  STATEMENT 

kk 

CAA 

kk 

CAA 

TI-ASC  ONLY 

:  USE  U-OPTION  WHEN  COMPILING  FOURIER  PACKAGES 

kk 

CAA 

kk 

CAA 

NX  IS  ARBITRARY  ( .GE.3) . 

kk 

CAA 

NY  MUST  BE 

OF  THE  FORM  PA2AAL,  WITH  L  *  1,  ...,  7. 

kk 

CAA 

USUALLY 

P*4,  5,  OR  &  - 

kk 

CAA 

HIGHER 

VALUES  OF  P  ARE  POSSIBLE,  BUT  ARE  LESS  EFFICIENT. 

kk 

CAA 

LQ  MUST  BE 

AT  LEAST  LQ1+LQ2  <SEE  CMNSUH). 

kk 

CAA 

kk 

CAA 

PROBLEMS  WITH  NON-ZERO  NEUMANN  OR  DIRICHLET  BOUNDARY  CONDITIONS 

kk 

CAA 

CAN  BE  SOLVED  BY  SUITABLY  MODIFYING  THE  BOUNDARY  SOURCE  TERMS. 

kk 

CAA 

< I.E.  MODIFY  S<I,J)  WHERE  1*1,  OR  I-NX,  OR  0*1,  OR  0*NY> 

kk 
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CAA  AA 
CAA  MUST  CALL  CMNSUH  FIRST  TO  SET  UP  CONSTANTS  REQUIRED  BY  SMNSUH.  AA 
CAA  IF  NONE  OF  ( NX f NY , DAX , DA Y , DAC  )  CHANGE,  THEN  SHNSUH  MAY  BE  AA 
CAA  CALLED  REPEATEDLY.  IF  ANY  OF  THESE  CHANGE,  CMNSUH  MUST  BE  AA 
CAA  CALLED  TO  RECALCULATE  THE  REQUIRED  CONSTANTS.  IF  TWO  OR  MORE  AA 
CAA  EQUATIONS  ARE  SOLVED  REPEATEDLY,  ONE  AFTER  THE  OTHER.  THEN  ONLY  AA 
CAA  ONE  CALL  TO  CMNSUH  FOR  EACH  IS  REQUIRED  IF  THEY  ARE  ALLOCATED  AA 
CAA  DISTINCT  CONSTANT  ARRAYS,  Q.  AA 
CAA  AA 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

C 

C  FORWARD  FFT. 

C 

CALL  COSODD ( S , W1 ,  U2 , Q , NDX , NX , NY ) 

C 

C  SOLVE  RESULTING  TRIDIAGONAL  SYSTEMS. 

C 

IQT* INT (Q( I ) )  +  1 
IQQ* IQT  +  3ANX 

CALL  TRANSF<W1 ,  NDX , NX , NY  ,  W2,NY) 

CALL  TGAUSS < W2 , W2 , NY ,  NY, NX,  Q <  IQT ) , Q <  IQQ  )  ) 

CALL  TRANSP < W2 ,  NY,  NY, NX,  Ml, NDX) 

C 

C  REVERSE  FFT. 

C 

DO  11  1*1, NX 

W1 ( I, 1 )*0.5AW1 (1,1) 

11  CONTINUE 

CALL  COSOD  I <  W1 , H ,  U2 , Q , NDX , NX , NY  ) 

RETURN 

C  END  OF  SMNSUH. 


r 


SUBROUTINE  CMMSSH ( DAX , DAY , DAC ,  NDX ,  NX ,  NY ,  QX , LQX , QY , LQY  ,  U) 

INTEGER  NDX , NX , NY ,  LQX, LQY 

DOUBLE  PRECISION  DAX(3,NX> ,DAY<3,NY> ,DAC 

DOUBLE  PRECISION  U 

DIMENSION  QX  <  NX , NX , 2 ) ,QY(LQY) ,  U(NX,1) 

C  DIMENSION  y ( NX , NY+5 ) 

C 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAA  AA 


CAA 

** 

CAA 

introduction: 

** 

CAA 

A* 

CAA 

COMPUTES  THE  CONSTANTS 

REQUIRED  BY  SMMSSH  TO  SOLVE  THE 

AA 

CAA 

FOLLOWING  (  MIXED-MIXED  >  HELMHOLTZ ' S  EQUATION  ON 

AA 

CAA 

AN  NXANY  MESH. 

AA 

CAA 

AA 

CAA 

FOR  I  *  1,  ...,  NX  AND  J  -  1,  ...,  NY 

AA 

CAA 

AA 

CAA 

AN<  J)  A 

H(I,  J*l>  + 

AA 

CAA 

AS (  J)  A 

H  < I ,  J-l)  + 

AA 

CAA 

AE ( I  >  A 

H( 1+1 , J  )  + 

AA 

CAA 

AW (I  )  A 

H  < 1-1 , J  )  ♦ 

AA 

CAA 

AP< I, J)  A 

H< I,  J  )  *  S< I,J> 

AA 

CAA 

AA 

CAA 

WHERE 

AA 

CAA 

AN<  J )  *  DAY  (  3  ,  J  ) 

J-l... NY-1, 

AA 

CAA 

0 

3*  NY. 

AA 

CAA 

AS ( J )  «  0 

J-  1. 

AA 

CAA 

DAY  <  1 ,  J ) 

J»2...NY, 

AA 

CAA 

AE ( I >  ■  DAX ( 3 , I > 

1*1. ..NX-1, 

AA 

CAA 

0 

I*  NX. 

AA 

CAA 

AW < I)  *  0 

I*  1. 

AA 

CAA 

DAX (1,1) 

1*2. ..NX, 

AA 

CAA 

AP< I, J>  ■  DAC+DAX < 2 

, 1 > - <  BN  <J)+BS(J>+BE( I)+BW< I) > 

AA 

CAA 

1*1. ..NX,  J«1...NY. 

AA 

CAA 

BN  <  J )  -  DAY (3,3) 

J-l. . .NY. 

AA 

CAA 

BS ( J)  -  DAY ( I , J ) 

J*l. . .NY. 

AA 

CAA 

BE ( I )  -  DAX ( 3 , I > 

1*1. ..NX. 

AA 

CAA 

BW< I)  -  DAX (1,1) 

I»l.. .NX. 

AA 

CAA 

AA 

CAA 

USUALLY  DAX ( 2 , I ) *0 . DO 

AND  DAY ( 2 , J) *0. DO . 

AA 

CAA 

AA 

CAA 

FOR  DIRICHLET-MIXED 

PROBLEMS  DAX(1,1>,DAX(3,NX>  ARE 

AA 

CAA 

SIMILAR  TO  OTHER  DAX  TERMS. 

AA 

CAA 

FOR  NEUMANN-MIXED 

PROBLEMS  DAX(1,1)*DAX(3,NX)*0.D0  . 

AA 

CAA 

AA 

CAA 

FOR  MIXED-DIRICHLET  PROBLEMS  DAY ( 1 , 1 > , DAY < 3 , NY >  ARE 

AA 

CAA 

SIMILAR  TO  OTHER  DAY  TERMS. 

AA 

CAA 

FOR  MIXED-NEUMANN 

PROBLEMS  DAY (1,1)*DAY(3,NY)*0.D0  . 

AA 

CAA 

AA 

CAA 

AA 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


CAA 

AA 

CAA 

ARGUMENT  LIST: 

SEE  OTHER  COMMENTS  FOR  DETAILS 

AA 

CAA 

AA 

CAA 

DAX,DAY,DAC  - 

HELMHOLTZ  COEFFICIENTS 

AA 

CAA 

AA 

CAA 

NDX 

1ST  ARRAY  DIMENSION  OF  S  AND  H  IN  THE  SOLVER 

AA 

CAA 

AA 

CAA 

NX, NY 

DIMENSIONS  OF  MESH 

AA 

CAA 

AA 

CAA 

QX  <  NX , NX , 2 )  - 

ARRAY  FOR  REQUIRED  CONSTANTS  (SEE  'USAGE'  BELOW) 

AA 

CAA 

AA 
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CAA  QY(LQY)  -  ARRAY  FOR  REQUIRED  CONSTANTS  AA 
CAA  AA 
CAA  LQX  -  SIZE  OF  OX,  MUST  BE  AT  LEAST  2ANXANX  AA 
CAA  AA 
CAA  LQY  -  SIZE  OF  QY ,  MUST  BE  AT  LEAST  NYA(NX+3>  AA 
CAA  AA 
CAA  W<NX,NY+3)  -  DOUBLE  PRECISION  UORKSPACE  ARRAY  AA 
CAA  AA 
CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAA A AAA AAAAAAAAAAAA 
CAA  AA 
CAA  SPECIAL  FEATURES  OR  NON-ANSI  USAGE  AA 


QX  IS  USED  AS  (DOUBLE  PRECISION)  UORKSPACE  IN  'ESYSON'. 
THEREFORE  ADVISABLE  TO  EQUIVALENCE  QX<1,1,1)  TO  A 
DOUBLE  PRECISION  VARIABLE. 


IT  IS 


CAA  AN  ERROR  MESSAGE  WILL  BE  OUTPUT  TO  CHANNEL  6,  AND  THE  PROGRAM  AA 
CAA  TERMINATED,  IF  ANY  OF  THE  FOLLOWING  CONDITIONS  ARE  DETECTED:  AA 
CAA  NX  LESS  THAN  3  AA 
CAA  NY  LESS  THAN  4  AA 
CAA  NY  LESS  THAN  NX  AA 
CAA  LQX, LQY  TOO  SMALL  AA 
CAA  DAX,DAY,DAC  NOT  CONSISTANT  WITH  HELMHOLTZ'S  EQUATION  AA 
CAA  AA 
CAA  THE  UORKSPACE  ARRAY,  U,  NEED  NOT  BE  DISTINCT  FROM  THE  AA 
CAA  UORKSPACE  SUPPLIED  TO  THE  SOLVER.  HOWEVER  THE  LATTER  MAY  AA 
CAA  CONSIST  OF  SINGLE  PRECISION  ARRAYS.  AA 
CAA  AA 
CAA  FOR  MORE  DETAILED  DOCUMENTATION  -  SEE  SMMSSH.  AA 
CAA  AA 
CAAAAAAAAAAAAAAAAAAAA A AAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAA 
CAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAA AAA AAA AAA AAA AAAAAAAAAAAAAAAAAAAAAAAAA 


END  OF  CMMSSH. 
END 


ni  <j. 


SUBROUTINE  SMMSSH< S ,H , NDX ,NX  ,NY  ,  QX , LQX , QY , LQY ,  Wl) 
INTEGER  NDX , NX , NY ,  LQX , LQY 
DIMENSION  S(NDX,NY) ,H(NDX,NY> 

DIMENSION  QX ( NX ,NX , 2 ) , QY  <  LQY  > 

DIMENSION  U1  ( NDX ,  NY ) 


Ckkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


Ckk 

kk 

Ckk 

DATE: 

DECEMBER  1980 

kk 

Ckk 

kk 

Ckk 

AUTHOR: 

A . J .  WALLCRAFT 

kk 

Ckk 

kk 

Ckk 

archivist: 

A . J .  WALLCRAFT 

kk 

Ckk 

kk 

Ckk 

ADDRESS: 

NORDA 

kk 

Ckk 

CODE  322 

kk 

Ckk 

NSTL  STATION,  MS  39529 

kk 

Ckk 

kk 

Ckk 

PHONE: 

601-688-4835  (COMMERCIAL) 

kk 

Ckk 

494-4835  (FTS) 

kk 

Ckk 

485-4835  (AUTOVON) 

kk 

Ckk 

kk 

Ckk 

DATE  OF  CURRENT  SEQUENCE:  DECEMBER  1982 

kk 

Ckk 

kk 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 

CAA 


INTRODUCTION: 

SOLVES  FINITE  DIFFERENCE  HELMHOLTZ'S  EQUATIONS: 

OVER  A  RECTANGLE, 

ON  A  STRETCHED- STRETCHED  MESH, 

WITH  EIGTHER  (ZERO)  NEUMANN  BOUNDARY  CONDITIONS  1/2  GRID 
DISTANCE  OUTSIDE  THE  MESH  OR  (ZERO)  DIRICHLET  BOUNDARY 
CONDITIONS  1  GRID  DISTANCE  OUTSIDE  THE  MESH  AT  J*1  AND  J*NY, 
AND  EIGTHER  (ZERO)  NEUMANN  BOUNDARY  CONDITIONS  1/2  GRID 
DISTANCE  OUTSIDE  THE  MESH  OR  (ZERO)  DIRICHLET  BOUNDARY 
CONDITIONS  1  GRID  DISTANCE  OUTSIDE  THE  MESH  AT  1*1  AND  I*NX . 

OR  ANY  EQUIVALENT  SET  OF  EQUATIONS;  FOR  EXAMPLE: 

HELMHOLTZ'S  EQUATION  IN  (SOME)  NON-CARTESIAN  COORDINATE 
SYSTEMS,  OR  MORE  GENERAL  (SEPARABLE)  ELLIPTIC  EQUATIONS. 

USING  THE  MATRIX  EIGENVALUE  DECOMPOSITION  METHOD. 

BOTH  THE  INPUT  ARRAY  S(NDX,NY>  AND  THE  OUTPUT  ARRAY  H ( NDX ,  NY ) 
ARE  ASSUMED  TO  CONTAIN  ACTIVE  SUBARRAYS  FROM  1  TO  NX  AND 
1  TO  NY.  THE  NON-ACTIVE  SUBARRAY,  NX+1 
OF  S  CAN  BE  ARBITRARY  ON  INPUT  AND  THAT 
INTACT. 


TO  NDX  AND 
OF  H  IS 


1  TO  NY, 
RETURNED 


IN  THE  CALLING  SEQUENCE  THE  INPUT  AND  OUTPUT  ARRAYS  MAY  BE 
SAME  ARRAY.  IF  THEY  ARE  NOT  THE  SAME,  THE  INPUT  ARRAY  IS 
RETURNED  INTACT  BY  THE  SOLVER. 


THE 


IT  IS  ASSUMED  THAT  NECESSARY  CONSTANTS  HAVE  BEEN  PRECOMPUTED 
EXTERNALLY  AND  ARE  STORED  IN  THE  ARRAYS  QX  AND  QY,  WHICH 
IMPLICITLY  DEFINES  THE  HELMHOLTZ  COEFFICIENTS  (DAX,  DAY,  DAC) 
TO  THE  SOLVER. 


FOR  1*1, 
SATISFY: 


NX  AND 


1, 


NY  THE  SOLUTION  WILL 


k  A 
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CAA 

AA 

CAA 

AN (  J  )  A 

H(I,  3*1)  + 

AA 

CAA 

AS (  J)  A 

H ( I ,  J-l)  + 

AA 

CAA 

AE ( I  )  A 

H  (  1  + 1 , J  )  + 

AA 

- - 

CAA 

AW ( I  )  A 

H(  1-1 , J  )  + 

AA 

* 

CAA 

AP ( I , J )  A 

H( I,  3  >  = 

S( I, J) 

AA 

CAA 

AA 

CAA 

WHERE 

A  A 

CAA 

AN ( J )  =  DAY ( 3 , J ) 

3=1.. .NY-1, 

AA 

. 

CAA 

0 

J*  NY. 

AA 

CAA 

AS ( J )  =0 

3*  1. 

AA 

CAA 

DAY ( 1 , J ) 

3=2. . .NY, 

AA 

CAA 

AE( I)  =  DAX ( 3 , I) 

II 

• 

• 

• 

7C 

X 

1 

H 

AA 

. ", 

CAA 

0 

I*  NX. 

A  A 

CAA 

AW( I)  >0 

1=  1. 

AA 

CAA 

DAX (1,1) 

1*2. . .NX, 

AA 

CAA 

AP ( I , J  )  =  DAC  +  DAX ( 

2, 1 )-(BN( 3>+BS( 3>+BE( I)+BW( I) ) 

AA 

CAA 

1=1. ..NX,  J= 1 . . . NY . 

AA 

CAA 

BN  (  J  >  =  DAY  ( 3 ,  J  > 

J=1 . . .NY 

AA 

CAA 

BS(J>  ■  DAY ( 1 , J ) 

3*1. . .NY 

AA 

CAA 

BE ( I)  =  DAX (3,1) 

1=1. ..NX. 

AA 

V 

CAA 

BW( I)  *  DAX (1,1) 

1=1. . .NX. 

AA 

CAA 

AA 

CAAAAAAAAAAAAAAAAA AAA AAA A AAA A AAA A A AAA A A A A A A A AAAAAAAAAAAAAAAAAAAA A A AAAAA A 

v 

CAA 

AA 

CAA 

LIBRARIES  REQUIRED: 

AA 

CAA 

• 

AA 

CAA 

LINPACK  BY 

AA 

CAA 

TRIDIAGONAL  LINEAR  EQUATION  SOLVER 

LIBRARY  BY  A.J.  WALLCRAFT 

AA 

CAA 

ROUTINE  HQR2  IS  FROM 

EISPACK 

AA 

CAA 

AA 

CAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAA A AAAAA  vaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 


ARGUMENT  LIST:  SEE  OTHER  COMMENTS  FOR  DETAILS 


S  <  NDX , NY ) 


INPUT  ARRAY 


H<NDX,NY>  -  OUTPUT  ARRAY 


NDX, NY 


-  SIZES 


AND  H 


-  ACTIVE  SUBDIMENSION 


AND  H 


QX ( NX  ,  NX  ,  2 )  -  CONSTANT  ARRAY,  PRECOMPUTED  EXTERNALLY 


QY  <  LQY ) 


-  CONSTANT  ARRAY,  PRECOMPUTED  EXTERNALLY 


-  DIMENSION  OF  QX 


-  DIMENSION  OF 


VII  < NDX , NY )  - 


WORKSPACE  ARRAY 
RETURNED  AND 


AVAILABLE  AFTER  EACH  CALL. 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


SPECIAL  FEATURES 


NON-ANSI  USAGE: 


DOUBLE  PRECISION  LIBRARIES  USE 


IMPLICIT  STATEMENT 


REPLACE  MATRIX-MATRIX  MULTIPLIES  (BELOW)  WITH  FASTEST  (MACHINE 
CODE)  VERSION  AVAILABLE  FOR  TARGET  COMPUTER.  THIS  MAY  INVOLVE 
USING  THE  TRANSPOSE  OF  QX  -  IF  SO  MODIFY  CMMSSH  APPROPRIATELY. 


CAA  AA 
CAA  NX  IS  ARBITRARY  (.GE.3)  BUT  MUST  NOT  BE  GREATER  THEN  NY.  AA 
CAA  NY  IS  ARBITRARY  (.GE.4).  AA 
CAA  LQX  MUST  BE  AT  LEAST  2ANXANX.  AA 
CAA  LQY  MUST  BE  AT  LEAST  NYA(NX+3).  AA 
CAA  AA 
CAA  PROBLEMS  WITH  NON-ZERO  NEUMANN  OR  DIRICHLET  BOUNDARY  CONDITIONS  AA 
CAA  CAN  BE  SOLVED  BY  SUITABLY  MODIFYING  THE  BOUNDARY  SOURCE  TERMS.  AA 
CAA  <I.E.  MODIFY  S(I,3)  WHERE  1=1,  OR  I=NX,  OR  3=1,  OR  3=NY>  AA 
CAA  AA 
CAA  MUST  CALL  CMMSSH  FIRST  TO  SET  UP  CONSTANTS  REQUIRED  BY  SMMSSH.  AA 
CAA  IF  NONE  OF  ( NX . NY , DAX , DAY , DAC >  CHANGE,  THEN  SMMSSH  MAY  BE  AA 
CAA  CALLED  REPEATEDLY.  IF  ANY  OF  THESE  CHANGE.  CMMSSH  MUST  BE  AA 
CAA  CALLED  TO  RECALCULATE  THE  REQUIRED  CONSTANTS.  IF  TWO  OR  MORE  AA 
CAA  EQUATIONS  ARE  SOLVED  REPEATEDLY,  ONE  AFTER  THE  OTHER,  THEN  ONLY  AA 
CAA  ONE  CALL  TO  CMMSSH  FOR  EACH  IS  REQUIRED  IF  THEY  ARE  ALLOCATED  AA 
CAA  DISTINCT  CONSTANT  ARRAYS,  QX,QY.  IF  TWO  OR  MORE  EQUATIONS  AA 
CAA  DIFFER  ONLY  IN  ( NY ,DAY , AC ) ,  AND  NOT  IN  <NX,AX>.  THEN  THEY  CAN  AA 
CAA  SHARE  A  SINGLE  QX  ARRAY,  ONLY  DISTINCT  QY  ARRAYS  ARE  REQUIRED.  AA 
CAA  AA 


CAAAAAAAAAAAA AAAAAAAAAAAAAAAAA A A AAA A A A A A A A A A AAAAAAAAAAAAAA AAA AAA A AAA AAA A 
CAAAAAAAAAAAAAAAAAAAAAAAA A A AAAAA AAAAAA A A A AAAAAAAAAAAAAAAAAAAAAAAA AAA A AAA 
C 

C  FORWARD  TRANSFORM  (BY  MATRIX-MATRIX  MULTIPLY). 

C 

DO  11  3=1, NY 
DO  11  1*1, NX 
W1(I,J)=0.0 
DO  11  K»1 , NX 

W1 ( I, J)*W1 (1,3)  ♦  QX(K, 1,2) AS (K, 3) 

11  CONTINUE 
C 

C  SOLVE  RESULTING  TRIDIAGONAL  SYSTEMS. 

C 

IQY=3ANY+1 

CALL  TGAUSS (W1,W1,NDX,NX,NY,  QY,QY ( IQY  >  ) 

C 

C  REVERSE  TRANSFORM  (BY  MATRIX-MATRIX  MULIIPLY). 

C 

DO  21  3*1, NY 
DO  21  1*1, NX 
H< I, 3) *0.0 
DO  21  K=1 , NX 

H( I,3)«H( 1,3)  ♦  QX(K, I,1)AW1(K,3> 

21  CONTINUE 
RETURN 

C  END  OF  SMMSSH. 

END 
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APPENDIX  B 


Example  of  Quasi -Geostrophic  Regional  Energetics  Graphical  Output 


GYRE  H0LLRND 


1  JULES  TWMSFER  UNITS  -  l.OE  05  J0ULES/3EC 


0.0  360.0  720.0  1060.01140.01800.02160.02520.02880.03210.03600.0 

TIME  IN  DAYS 


XX). 02160.02520. 02880.03240. 03600.0 
IN  DAYS 
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